An electronic literature search was performed in Medline and Embase databases to identify all eligible studies on the association between pre-diagnostic 25(OH)D concentrations and RA risk that were published from inception to 21 December 2021. The search strategy included keywords and MeSH terms covering a comprehensive list of possible variations of ‘vitamin D’ and its forms, ‘rheumatoid arthritis’, and prospective study designs. The full search strategy for Medline and Embase is listed in Additional File 1. The systematic review protocol was registered with Prospero on 2 August 2021 (CRD42021262855) [20].

Results from the database searches were merged and duplicates removed using Covidence, a web-based software platform for systematic review management [21]. Records were independently screened by two investigators (JLC and RC) based on titles and abstracts in the first step. The full texts of the remaining records were then independently screened for eligibility in duplicate and uncertainties regarding inclusion were resolved by discussion with all investigators. The reference lists of included papers were manually searched and a forward citation search was performed to ensure no potentially relevant studies were missed. The inclusion criteria were as follows: prospective cohort studies, case-cohort studies, nested case–control studies, or case–control studies with retrospective collection of exposure data which reported adjusted relative risk (RR) estimates (RRs, hazard ratios (HRs), odds ratios (ORs), incidence rate ratios) and 95% confidence intervals (CIs) for RA risk in relation to 25(OH)D; with 25(OH)D concentrations measured in blood samples collected prior to diagnosis of RA; and participants aged 18 years and older. The outcome of interest was RA diagnosis, by self-report or clinician report, using American College of Rheumatology (ACR), European Alliance of Associations for Rheumatology (formerly European League Against Rheumatism, EULAR) [22], or International Classification of Diseases (ICD) criteria. Conference abstracts were excluded because they are unlikely to contain all the information required for statistical analyses and for evaluation of study quality. If there were multiple publications from the same study population, only the most recent was included. Reviews, letters, comments, editorials, meta-analyses, meta-syntheses, prediction models, protocols, cross-sectional studies, case reports, and epidemiological studies in patient populations were excluded. We also manually excluded non-human studies and non-English publications.

Data extracted from the eligible papers included: first author name, year of publication, country/location, study name or description, study period, duration of follow-up, number of participants or controls, number of cases, age and sex distribution, study design, outcome assessment, method of 25(OH)D measurement, mean/median 25(OH)D, comparison (the contrast or metric of 25(OH)D), RR/OR/HR and 95% CI, and any adjustments made for covariates. If several multivariable models were reported, the most fully adjusted effect estimate was extracted and used for the meta-analysis. If an article reported pooled results from more than one study, only the individual results were used in this meta-analysis. Data were extracted by one investigator (RC) and separately verified by two investigators (JLC and AKH).

The quality of each of the studies included in the meta-analysis was assessed with the Newcastle–Ottawa scale [23]. Each study was appraised using the categories of selection of cases and controls, comparability, and ascertainment of the exposure for case–control studies, and participant selection, comparability, and ascertainment of the outcome of interest for cohort and nested studies, with a maximum possible score of 9 stars. For publications reporting on multiple studies, each study was rated separately.

Prior to meta-analysis, results from all studies were converted to estimates for a 25 nmol/L increment in 25(OH)D. If a study reported results for 25(OH)D in ng/ml, concentrations were first converted to nmol/L by multiplying by 2.496. Given the low prevalence of RA, different effect measures (HRs, ORs) were considered equivalent to RRs and used directly. If a study reported results separately by sex, or other subgroups, but not overall, the subgroup-specific estimates were pooled using a fixed-effect model before inclusion in the meta-analysis. If a study reported a CI with insufficient precision (due to rounding), the CI was estimated using the reported p-value. Because few studies met the eligibility criteria for inclusion, we used a Bayesian meta-analysis approach to estimate the summary RR and 95% credible interval (CrI) for RA in relation to a 25 nmol/L increment in 25(OH)D concentration. The CrI, similar to a frequentist CI, quantifies the uncertainty around the point estimate. Heterogeneity of RRs across studies is to be expected, given the differences in populations, study designs, and methods for ascertainment of outcomes. A frequentist random-effects model can reliably handle heterogeneity in a meta-analysis of a large number of studies, however, these models are not well-suited to meta-analyses of few studies. A Bayesian approach to meta-analysis accounts for uncertainty in the estimate of heterogeneity inherent in such situations [24]. Additionally, a prior distribution for τ can be used to stabilize the heterogeneity estimate around a reasonable value. Empirically derived priors for τ, based on previously published meta-analyses, have been developed, and here we employed the prior based on studies examining non-pharmacological exposures and onset of a new chronic disease [25]. We reported I2, the percent of total variation that is due to true between-study heterogeneity, which is derived from the posterior distribution of τ. The prior distribution for the log RR is weakly informative (normal with mean 0 and standard deviation 0.56), based on the assumption that the RR will have a magnitude of less than 3.

To investigate potential heterogeneity in the results, subgroup analyses were conducted by sex among studies that reported sex-specific estimates. To evaluate whether the results were driven by one large study or a study with an extreme result, a sensitivity analysis was performed sequentially excluding each study from the analysis one at a time and assessing the impact on the summary estimate. As a further sensitivity analysis, frequentist random-effects models were run to estimate the summary RR and 95% CI for risk of RA per 25 nmol/L increment in 25(OH)D. Risk of publication bias was evaluated by visually inspecting funnel plots and using the Egger linear regression test [26] and Begg rank correlation test [27].

All statistical analyses were performed using R version 4.1.0 [28], and the metafor version 3.0–2 [29] and bayesmeta version 2.7 [30] packages.